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Abstract. Multi-mode entanglement is investigated in the system composed of N 
coupled identical harmonic oscillators interacting with a common environment. We 
treat the problem very general by working with the Hamiltonian without the rotating- 
wave approximation and by considering the environment as a non-Markovian reservoir 
to the oscillators. We invoke an iV-mode unitary transformation of the position and 
momentum operators and find that in the transformed basis the system is represented 
by a set of independent harmonic oscillators with only one of them coupled to the 
environment. Working in the Wigner representation of the density operator, we find 
that the covariance matrix has a block diagonal form that it can be expressed in terms 
of multiples of 3 x 3 and 4x4 matrices. This simple property allows to treat the 
problem to some extend analytically. We illustrate the advantage of working in the 
transformed basis on a simple example of three harmonic oscillators and find that the 
entanglement can persists for long times due to presence of constants of motion for 
the covariance matrix elements. We find that, in contrast to what one could expect, 
a strong damping of the oscillators leads to a better stationary entanglement than in 
the case of a weak damping. 
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1. Introduction 

Controlled dynamics and preservation of an initial entanglement encoded into a 
continuous variable system of harmonic oscillators coupled to a noisy environment are 
challenging problems in quantum information technologies [U |2] . The coupling induces 
decoherence phenomena, such as decay and dissipation that reduce and even can destroy 
the initial entanglement over a finite evolution time [3j HJ IS]- Dynamics of an open 
quantum system are usually studied in terms of the master equation of the reduced 
density operator whose structure depends on the nature of the environment to which 
the system is coupled. It has been noted that the dynamics crucially depend on whether 
the oscillators interact with a common or independent local environments. In the later 
case the interaction usually leads to a degradation of the entanglement whereas in the 
former, the environment can not only create decoherence, as it usually does, but may 
act as a source of coherence that not only preserves the initial entanglement but also 
creates an additional entanglement. A series of papers accounts these properties for the 
case of two coupled harmonic oscillators being in contact with a Markovian thermal 
reservoir and the work of Liu and Goan [6], Maniscalco et al. [7] and Horhammer 
and Biittner [8] accounts for a non-Markovian thermal bosonic reservoirs. Detailed 
discussions and extensive reference lists devoted to the decoherence of two harmonic 
oscillators can be found in Refs. P HQl HH H21 H3l HH US! UHl HTl UHl [I9l [20] EI]. Non- 
Markovian quantum dynamics of open systems has been discussed by others, notably by 
Breuer and Vacchini [22] , who provide the memory kernel treatment and illustrate it for 
various examples and applications. In all these studies a general conclusion made is that 
entanglement dynamics depends on the form of the reservoir and the non-Markovian 
nature of the reservoir preserves entanglement over a longer time. 

More important in the quantum technologies is the characterization and the study 
of dynamics of a large number of harmonic oscillators that are crucial for the study of 
quantum coherence, entanglement, fluctuations and dissipation of mesoscopic systems. 
The correct understanding of the mechanism responsible for entanglement evolution in 
the system is essential for designing iV-atom systems for quantum information processing 
and quantum computation. The key problem is to find the master equation for N 
harmonic oscillators coupled to an environment that can be solved in a simple and 
effective way. How to treat such a composed system in the most effective way and how 
to understand its complicated dynamics are challenging questions that still have not 
been resolved. 

In this paper, we pursue a research that especially addresses these questions. In the 
approach, we treat the problem very general, fully accounting the non-RWA dynamics 
and considering the environment as a non-Markovian reservoir to the oscillators. We 
introduce an iV-mode unitary transformation of the position and momentum operators 
and find that in the transformed basis the system is represented by a set of independent 
harmonic oscillators with only one of them coupled to the environment. This fact makes 
the problem remarkably simple that the relaxation properties of N harmonic oscillators 
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follow the same pattern as a single harmonic oscillator. This property also leads to two 
different time scales of the evolution of the system, a short time scale where the dynamics 
are strongly affected by the relaxation process, and a free of the relaxation long time 
scale. We distinguish those two time scales by working within the correlation matrix 
representation, also known as the covariance matrix. We also consider squeezing of the 
position and momentum variances as a practical measure of a three-mode entanglement. 
We compare squeezing with the negativity [231 121] and show that suitably transformed 
(rotated) quadrature components of the field modes exhibit squeezing whenever there 
is entanglement between the modes and vice versa. 



2. The model 



We consider a system composed of N mutually coupled identical harmonic oscillators of 
mass M and frequency Q that are simultaneously interacting with a common thermal 
bath environment (reservoir). The system is determined by the Hamiltonian, which in 
terms of the position qi and momentum pi operators can be written as 



H = H S + H £ + V S + V, (1) 



where 



N ' 2 1 
Pi , 1 



X 2M 2 

is the free Hamiltonian of the harmonic oscillators, 

is the Hamiltonian of the common reservoir to which the oscillators are coupled, 

N 

^ = A EE^- ( 4 ) 

i=l j>i 

is the interaction between the oscillators, and 

N 

V = x ^q n qi (5) 

n i=l 

is the interaction between the oscillators and the reservoir. 

In equations ([I])-©, the parameter A stands for the coupling constant between the 
oscillators, and A n is the coupling strength of the oscillators to the reservoir. We model 
the environment as an ensemble of harmonic oscillators of mass m n and frequency u n 
that interact bilinearly through their position operators q n with the oscillators. 

The system of harmonic oscillators coupled to an environment is usually described 
in terms of a reduced density operator p, which is obtained by tracing the density 
operator of the total system over the reservoir operators. Instead of working in the bare 
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basis, (qi,pi 
operators 



we introduce an iV-mode unitary transformation of the systems' position 



Qk 



N-k 
N-k+1 

N 



k=l 



Qk 



N-k 



N 

i=k+l J 



fc = l,2,3,...,JV-l, 



(6) 



and the same for the momentum operators. We note that the transformations involve 
anti-symmetrical (qi,pi) and symmetrical (c[n,Pn) combinations of the position and the 
momentum operators, a close analog of the symmetric and antisymmetric multi-atom 
Dicke states [251 HH 123 • 

In order to derive the master equation for the density operator p of the system, 
we use the standard method involving the Born approximation that corresponds to 
the second-order perturbative approach to the interaction between the oscillators 
and the environment, but we do not make the rotating- wave (RWA) and Markovian 
approximations. We find that in terms of the transformed operators the reduced density 
operator p satisfies the master equation 



i 

h 



H s + -Mn 2 N (t)q 2 N ,P 



jlif) [Qn,{Pn,p}} 



D(t)[q N , [qN,p\] - t/(*)[«at, [Pn,p]\ 



in which the Hamiltonian H s of the coupled oscillators is of the form 



N 



2M 2 lHl 



(7) 



(8) 



where 



M 



1,2, 



iV-1, 



n 



N 



n> + <w-i)£ 



(9) 



are the effective frequencies of the oscillators. Note that the frequency Qn of the 
oscillator coupled to the environment differs from that of the remaining independent 
oscillators. It means that the reservoir affects the evolution of only one of the oscillators 
leaving the remaining N — 1 oscillators to evolve freely in time. The dynamics of the 
N— th oscillator that is affected by the reservoir are determined by the following time- 
dependent coefficients 

2 



M 



dt\ cos(Ojv*i)n(*i 



(10) 



represents a shift of the frequency of the oscillator due to the interaction with the 
environment. It includes the frequency renormalization that leads to a finite Lamb 
shift [21. 
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The time dependent parameter 

i r* 

lN(t) = —— dhsmiQ^Uih) (11) 



is the dissipation coefficient, and 



1 f l 

D N (t) = - / d*icos(n^i)K*i), ( 12 ) 
1 /•* 

■M*) =-1777- / *isin(^i)K*i), (13) 

are diffusion coefficients. 

The time dependent functions H(t) and u(t) appear as the dissipation and noise 
kernels, respectively, and are given by 

1 C°° 
U ( t ) = ^Y t ^n(t),q n (0)]) = J q du J (u) Sm(u t) , (14) 

1 f°° 
=^Y, X l({ ( ln(t),qn(0)}} = J Q duJ(u)cos(ujt)[l + 2N(uj)i (15) 

where J{oj) is the spectral density of the modes of the environment. For a Gaussian-type 
spectral density 

■7M = Vg^e"^, (16) 

where A is cut-off frequency that represents the highest frequency in the environment, 70 
is proportional to the coupling strength between the N— th oscillator and the 
environment, and n determines the type of the reservoir. For n — 1, the environment 
is called an Ohmic reservoir, for n > 1 it is called supra-Ohmic, whereas for n < 1 it is 
called a sub-Ohmic reservoir. 

It should be stressed here that the derivation of the master equation holds under 
the Born approximation which takes the interaction between the oscillators and the 
reservoir only to the second order of the coupling strength A n . 

In the transformed basis, the Hamiltonian of the system exhibits interesting 
properties. First of all, we observe that the system is represented by a set of N 
independent oscillators with only one of them being coupled to the environment. The 
oscillator effectively coupled to the environment is that one corresponding to the 
symmetric combination of the position and momentum operators. In addition, the 
effective frequency Qn of the oscillator coupled to the environment differs from that of 
the remaining independent oscillators. The oscillators effectively decoupled from the 
environment may be regarded as composing a relaxation-free subspace. It should be 
stressed that the subspace is not a decoherence-free subspace. We shall demonstrate 
that the subsystem of the "relaxation free" oscillators still can evolve in time that may 
lead to decoherence. We will recognize that only a part of the subspace can be regarded 
as a decoherence-free subspace. 
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3. Covariance matrix 

We study dynamics of the system in terms of the Wigner characteristic function, which 
for an iV-mode Gaussian state can be written in terms of a covariance matrix as [29] 

X (X) = exp f-ljCVxA , (17) 

where X = col(gi,pi, q2,p2, ■ ■ ■ , Qn,Pn) is an 2N dimensional column vector of the 
transformed operators, and V is the covariance matrix whose elements are defined as 



V iJ = Ti({AX l ,AX 3 }p), (It 



with 



{AX h AXj} = - (AXiAXj + AXjAXi) , (19) 



1 

2 

and Xi is the ith component of the vector X. 

The covariance matrix is composed of 4N 2 elements. However, due to the 
symmetrical property that V^- = Vji, it is enough to find the diagonal elements and 
those off-diagonal elements with i < j to completely determine the matrix. Thus, the 
number of elements that have to be found is equal to N(2N+1). Technically, it is done by 
using the definition (fl8l) and the master equation ([7]) form which one finds the equations 
of motion for the covariance matrix elements that then are solved for arbitrary initial 
conditions. However, the equations form a set of coupled linear differential equations 
whose number is large even for a small number of oscillators. Therefore, the dynamics of 
coupled harmonic oscillators have usually been studied by employing numerical methods. 

We propose a different approach that illustrates the advantage of working in the 
basis of the transformed position and momentum operators. As we shall see, the 
approach allows to determine the covariance matrix elements in an effectively easy way 
requiring to solve separate sets of equations composed of only a small number of coupled 
differential equations. 

From equation (fT81 and the master equation ([7]), we find a set of inhomogeneous 
differential equations for the covariance matrix elements, which can be written in a 
matrix form as 

V N (t) = c N (t)v N (t) + hP N (t), (20) 

where 

V]^(t) = Col(^n, V12, V22, • • • , ^2AT-l,2Ar-l) V2N-l,2Ni V2N,2N) (21) 

is a column vector of the covariance matrix elements, 

F N (t) = col(0, 0, 0, ... , -f N {t), 2hD N {t)) (22) 

is a column vector composed of the inhomogeneous time-dependent terms, and Cat(£) 
is an N(2N + 1) x N(2N + 1) block diagonal matrix of the time-dependent coefficients. 
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The matrix Cpf(t) is a direct sum of small size matrices 



C N (t) 



N 



n=2 



(A 1 (0) ® A[ N ~ n \0) ® A 3 (t) 



Mt), 



(23) 



where 



/ 



Ai(0) 



2M" 



—MVL 2 F 



\ 



o 



-2Mtt 2 F 



2M" 



(24) 



-MQ%{t) - 7 (t) 



\ 

M" 1 



(25) 



A 3 (0 



and 



A 4 (0) 



\ 





-2ikm^(t) -2 7 (t) y 




/ 





M" 1 M" 1 


o \ 




-MQ 2 N (t) 


- 7 (t) M- 1 




-MVtp 


M" 1 


V 





-MVL 2 F -MVL 2 N {t) - 


"7(0 / 


/ 





M- 1 M- 1 N 






-MQ 2 F 


M" 1 






—Mflp 


M- 1 




V 










(26) 



(27) 



with n 2 N (t) =Vl 2 N + Q 2 N (t) and 7 (t) = 2 7jv (t). The superscript (JV - n) in Aj w n) (0) is 
understood as the number of the A 4 (0) matrices appearing in the direct sum. Thus, for 
iV = 2, no matrix A 4 (0) is involved in Cj^(t), one matrix A 4 (0) is involved for N = 3, 
and so on. 

There are several interesting and important conclusions arising from equation ( 1231) . 
Firstly, the equations of motion group into decoupled subsets of smaller sizes involving 
only three and four equations. In other words, the block diagonal matrix Cat(^) is 
composed of 3 x 3 and 4x4 matrices. Secondly, the evolution of an iV > 2 system 
of harmonic oscillators is determined by the same matrices as that determining the 
evolution of N = 2 oscillators [9l HQl HH H2J UHl HH HS1 [III HZl HE! E0] - Thirdly, the 
matrices A^O) and A 4 (0) are independent of time. This means that the time evolution 
of the covariance matrix elements whose dynamics are determined by Ai(0) and A 4 (0) 
can be found in an exact analytical form. Fourthly, the matrices A^O) and A 4 (0) are 
independent of the relaxation coefficient 7 . Thus, they reflect features of the N — 1 
oscillators that are effectively decoupled from the environment, and as such could be 
regarded as determining a relaxation-free subspace. Finally, the matrices A 2 (i) and 
As(t) are explicitly dependent on time through the relaxation terms 7 (t). Therefore, 
they represent dynamics of the oscillator effectively coupled to the environment. The 
explicit time dependence of the matrices (1251) and ( 1261) results from the non-Markovian 
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nature of the environment, and the matrix becomes time independent in the case of a 
Markovian situation. Thus, in the case for a Markovian reservoir, all the covariance 
matrix elements can be found analytically. 

It should be noted here that the relaxation-free subspace cannot be regarded as 
a decoherence-free subspace because the covariance matrix elements can undergo a 
periodic time evolution that may result in a periodic decoherence. To explore this, 
we look into the properties of the matrix Ai(0). It is easy to note that the determinant 
of the matrix Ai(0) is equal to zero. Mathematically, it means that among the three 
matrix elements involved, Vn, V\ 2 and V 22l there is a linear combination whose equation 
of motion is decoupled from the remaining equations. It is easy to show that the linear 
combination 

V+ = MiljA]. + jjV 22 , 

V n = MQ 2 F V n - -^V 22 (28) 
obeys V£ = and the remaining elements form a set of two coupled equations 

v u =4n 2 F v 12 , V 12 = -V n , (29) 

where = MQ%V n - {l/M)V 22 . 

The property of V{i = indicates that the linear combination is a constant 
of motion, i.e. V^(t) = V^(0). In other words, V{[(t) does not change in time and 
retains its initial value for all times. Physically, if initially the system was prepared 
in a state such that V^O) 7^ and with the other elements of the covariance matrix 
equal to zero, it would remain in that state for all times. For example, if the initial 
state is an entangled state, the initial entanglement of the system will remain constant 
in time. Therefore, the subspace composed of the V{i(t) element can be regarded as a 
decoherence-free subspace. 

The remaining matrix elements V{[ and V\ 2 can undergo a temporal evolution. 
Since there is no damping involved in the equations of motion ( 1291) . the solution would 
lead to the matrix elements continuously oscillating in time. It is easy to find that the 
solution of equation (129]) has a simple form 

V n (t) = V n (0) cos(2fi F t) + 2Q F V 12 (0) sin(2fi F t), 

V 12 (t) = V 12 (0) cos(2ft F t) - %21 Bin(2fi F t), (30) 

from which we see the matrix elements continuously oscillate in time with frequency 2Qp. 
This indicates that the system will never reach a stationary time-independent state 
unless V{l(0) = Vi 2 (0) = 0. We stress that the continuous in time oscillations are not 
related to the non-Markovian nature of the reservoir as the matrix Ai(0) determines 
dynamics of the oscillators that are not coupled to the reservoir. 

It is also found that the determinant of the matrix A^O) is equal to zero. Thus, 
following the above analysis, we can find that the set of the equations of motion for 
the covariance matrix elements determined by the matrix Azt(0) can be reduced to 
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two constants of motion and two equations of motion with the same coefficients as 
in equation (1291) . 

On the basis of the above analysis, we may draw a conclusion that the set of 
equations of motion for the covariance matrix elements can be converted into (N — l) 2 
constants of motion and a smaller size set of coupled equations determined by a matrix 

C' N (t) = Al 2 j (0)©Af^(t)©A 2 (t), (31) 

where A 5 (0) is a 2 x 2 matrix composed of the coefficients of the two coupled equations 
of motion fl29|). 



4. Multi-mode entanglement and squeezing 

We have already shown that due to the presence of the constants of motion in the 
evolution of the covariance matrix elements, the dynamics of the system, even after 
a long time, may strongly depend on the initial state. Since we are interested in the 
evolution of an initial entangled state and it is well known that multi-mode squeezed 
states are examples of entangled states, We have already shown that due to the presence 
of the constants of motion in the evolution of the covariance matrix elements, the 
dynamics of the system, even after a long time, may strongly depend on the initial 
state. Since we are interested in the evolution of an initial entangled state and it is well 
known that multi-mode squeezed states are examples of entangled states, we consider 
two experimentally realizable initial squeezed vacuum states with markedly different 
squeezing behaviors. We also demonstrate that with the two specific initial states, the 
problem of treating the dynamics of N harmonic oscillators simplifies to analysis of the 
properties of only those constants of motion and the matrices which involve only the 
diagonal elements of the covariance matrix. 

In the first example, we consider the most familiar multi-mode continuous variable 
Greenberger-Horne-Zeilinger (GHZ) entangled state [301 EI] 



with 



U\ = exp 



(33) 



E( 4 ^H 6 ') 2 )- Hx - 

where r is the squeezing parameter and the ket |0f, 4 ) represents the state with zero 
photons in each of the N modes. This GHZ state for N = 3 has been realized 
experimentally by two groups [331 EI] • 

In the second example, we assume that the system is initially prepared in a pure 
non-symmetric multipartite squeezed state of the form 

Af 

i^H^ni ^ ( 34 ) 
i=i 
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where the squeezing transformation is of the form 



{N-l N-l } 

i=i i^j=i J 



(35) 



Here, each of the N — 1 relaxation-free modes is correlated to a degree ro with the 
damped mode, whereas the relaxation-free modes are correlated between themselves 
to a degree r s . Practical schemes for generation of such a state have recently been 
discussed [33]. For example, it could be created by use of concurrent interactions in 
a second-order nonlinear medium placed inside an optical resonator, which might be 
realized experimentally in periodically poled KTiOP04 [36J. We will use this example 
to demonstrate the dependence of stationary entanglement on the amount of correlations 
initially encoded into the relaxation-free modes. 

One manifestation of the squeezed properties of the states is entanglement between 
different modes. We examine this entanglement property shortly, but first we examine 
the manifestation of the squeezed correlations in the form of the covariance matrix. 

The choice of the initial state ( )32|) is a consequence of the diagonal form of the initial 
covariance matrix. In particular, the initial values of the covariance matrix elements of 
the two-mode case are 



Mo) 



VU(0) = ie 



-2r 



MO) 



Mo) = ^, 



(36) 



whereas for the three-mode case the initial elements are 

1 



Vix(0) 

Mo) 



V 33 (0) = V 66 (0) 

Mo) 



-2r 



MO) = 2^, 



With the asymmetric squeezed state (1341) . 
diagonal and has the following symmetric form 



(37) 

the initial covariance matrix is not 



V(0) 



/ Mo) 



Mo) 



Mo) 
\ o 





Mo) 



Mo) 



Mo) 



Mo) 
o 

Mo) 
o 

V 35 (0) 






Mo) 



Mo) 



Mo) 



Mo) 
o 

Mo) 
o 

Mo) 
o 



o \ 

Mo) 
o 

Mo) 
o 

Mo) J 



(3* 



where the explicit expressions for the non-zero matrix elements are given in the 
Appendix A. 

Before proceeding further with the analysis of the entangled and squeezing 
properties of the system, we return for a moment to the solutions for the covariance 
matrix elements. We have seen that with the states (132]) . the initial covariance matrix 
is diagonal. An immediate consequence of the diagonal form of the initial covariance 
matrix is that for t > 0, the diagonal elements will be different from zero and only those 
off-diagonal elements whose the equations of motion were coupled to the equations of 
motion for the diagonal elements. It is easy to see from equation f l20|) that the non-zero 
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elements are determined by the matrices Ax(0) and A 2 (t). Thus, dynamics of a system 
composed of N harmonic oscillators can be readily determined from properties of the 
two simple 3x3 matrices. 

For a multi-mode quantum state, if the variance of the quadrature operator g$ 
meets the inequality A(^) < 1/2, then we say the state exhibits ordinary multi-mode 
squeezing. The minimum variance corresponds to the optimal multi-mode squeezing. 
However, the ordinary squeezing is produced by both one and two-mode correlations, 
whereas entanglement is solely related to the two- mode correlations [32]. Thus, the 
ordinary squeezing does not necessarily mean entanglement. We may distinguish 
between the contributions of the one and two-mode correlations to the variances and then 
determine the multi-mode squeezing by performing suitable unitary transformations of 
the mode operators. 

We illustrate this procedure for the case of three modes since the GHZ state (1521) 
for N = 3 is an example of multipartite entangled state whose entanglement is 
shared by more than two parties. Moreover, the three-mode GHZ state has been 
realized experimentally [33], [31] and also has been successfully applied to demonstrate 
quantum teleportation [3U [37] and quantum dense coding |33j. We will demonstrate 
the equivalence between the three-mode squeezing and the negativity criterion for 
entanglement. The two-mode case, N = 2, has been extensively studied in the 
literature [3H]. First, we make a local squeezing transformation on each of the modes, 
which results in transformed annihilation operators of the form [39J 

a x = aie tf (uiU2 - e 2 ^" e) v x v 2 ) + a{e^ e { Ul v 2 - e~ 2i{ip - e) Vl u 2 ) , 

a 2 = a 3 = a 2 e ie { Ul u 2 + e 2 ^" e V 2 ) - a^V^ + e^^W), (39) 

where Ui = cosh(rj) and Vi = sinh(rj) (i = 1,2) and the transformation has been made 
with the squeezing parameter r% and the phase angle ip on the mode 1, and with r 2 and 
the phase angle 9 on the modes 2 and 3. 

We use the Wigner characteristic function, which in terms of the above specifically 
chosen transformation can be written in a Gaussian form as 

t) = exp 1 -J/IG/I 1, 1 , (40) 



2 

where \x = (yi,xi,y 2 ,x 2 ,y3,xs) is a vector composed of the real ijj and imaginary 
Xj (j = 1,2,3) parts of the phase-space variables corresponding to operator a-,-, and G 
is the correlation matrix of the form 

/a0c0c0\ 
b d d 



G 



c a c 
d b d 
c c a 
\0 d d b J 
Note, the matrix G involves only four parameters that are 
o=-2/ 3 e^, b = -2f 3 e~ 2r ^ 



(41) 
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where hi = (\fi 



fi 
h 
h 



2h 1 e 2r2 , d = 2h 2 e~ 2r \ (42) 
- f 2 ) and h 2 = ~(\fi\ + A), with 
= m 2 u\ + m* 2 e Aiv v 2 x + 2m 4 e 2l ' p u 1 v 1 , 
-- (m 2 e" 2 ^ + m;e 2 ^) u x v x + m 4 (l + 2vf), 

-- A\m 1 \u 1 v 1 +m 3 (l + 2vl), (43) 

and mj are linear combinations of the covariance matrix elements given in the bare 
basis 

mi = l -(V{ 1 -Vl 2 -2iV{ 2 ), 
^2 = ^3-^4-^(^4 + ^3), 

m 3 = V[ x + V£ 2 , m 4 = V( 3 + V 2 \. (44) 

The entangled nature of the three-mode squeezed states is clearly exhibited by the 
presence of the off-diagonal terms in the correlation matrix G. 

The squeezing parameters ri,r 2 and the phase angles <p, 9 appearing in the 
transformation of the field operators can be carefully chosen to match the form of the 
correlation matrix G with the form of the covariance matrix V in the bare basis. In this 
way we can achieve the equivalence between three- mode squeezing and entanglement. 
This can be done by choosing the squeezing parameters as 



,2ri 



with mi 



m 3 — 2\mi\ 
m 3 + 2|mi| 
\rrii\ exp(2i(p) and f\ = 



|fe 2 |+/3 

\hi\+f 3 



(45) 



|A|exp(2i0). 

Having available the time dependent solutions for the covariance matrix elements, 
we then can easily find the characteristic function that allows us to compute variances 
of the position operators and momentum operators 

3 



for k — 1, 2, and 
X 3 



3-k 
2(4 - k) 



3-k 



j=k+l 



-l\ 



3-k 
2(4 -k) 



E 



k ^ 3 

j=k+i 



H.c. 



1^ 



st 



Y 3 



j'=i 



(46) 



(47) 



The variances are involved in the criterion for multi-mode squeezing that fluctuations of 
the correlations between three modes are squeezed if and only if the sum of the variances 
((AXj) 2 ) and ((AY,) 2 ) with i ^ j satisfies the following inequality [23] 

((AX,) 2 ) + ((AY}) 2 > < 1, i,j = 1,2,3. (48) 

Among the permutations of the variances involved on the left-hand side of equation (l48p . 
there might be more than one satisfying inequality condition for multi-mode squeezing. 
In this case, we choose the combination that reflects the largest squeezing. 
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To quantify entanglement between the modes, we adopt the negative partial 
transpose criterion that is known as the necessary and sufficient condition for 
entanglement of two- and three-mode Gaussian states. We will compare the criterion 
with the squeezing criterion to quantify squeezing as an alternative necessary and 
sufficient condition for entanglement. The advantage of the squeezing criterion over 
the negativity is that the former can be directly measured in experiments whereas the 
later can be inferred from the reconstruction of the density matrix of the system. 

The partial transpose criteria are based on the non-positive partial transpose of a 
matrix [32] 

Wit^ + ^ia, j = 1,2,3,... (49) 

where Tj is the partial transpose matrix with the transposition made on the jth 
mode block and a is a block diagonal symplectic matrix. It has been shown that 
multi-mode Gaussian states are not completely separated when for all j there are 
negative eigenvalues of the matrix ()49l) . The eigenvalues can be degenerated or non- 
degenerated. However, for a system of identical oscillators the covariance matrix V'(t) 
is permutational symmetric, so all the negative eigenvalues are degenerated. We denote 
them by a parameter 7/_ and call it as the negativity criterion for entanglement. 

5. Temporal evolution of squeezing and entanglement 

We now perform numerical analysis of time evolution of multi-mode squeezing and 
entanglement in a system of two and three mutually interacting harmonic oscillators 
simultaneously coupled to an environment. We will illustrate the advantage of working 
in the transformed basis to obtain a simple interpretation of the results. In particular, 
to understand short time non-Markovian dynamics of entanglement and to provide 
conditions for optimal and stable long time entanglement. In addition, we compare 
the time evolutions of the variances and the negativity to find if the condition for 
three-mode squeezing could be used as the necessary and sufficient condition for three 
mode entanglement. In all cases considered here, we assume that the oscillators interact 
with an Ohmic reservoir (n = 1) of temperature fc#T = 10HQ with the Boltzmann 
distribution of photons characterized by the mean occupation number N(Q) = 9.5083. 

We first consider the case of mutually independent oscillators with A = 0, but 
interacting with the environment. Figure [1] shows the negativity and variances as a 
function of time for the initial symmetric squeezed state with different degree 
of squeezing r. First of all, we note that at times where squeezing occurs there is 
entanglement, and vice versa, at times where entanglement occurs, there is squeezing. In 
addition, we see a threshold value for the degree of squeezing r at which a continuous in 
time entanglement occurs. The threshold that corresponds to entanglement undergoing 
the phenomenon of sudden death, occurs at r = 1.498. It is interesting to note that the 
same threshold value for r has been predicted for the two-mode case [T7] . 

The presence of the threshold value for r at which continuous in time entanglement 
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Figure 1. Time evolution of the negativity -q^ and the combined variance ((AXi) 2 ) + 
((AF 3 ) 2 ) for 7 = 0.05, A = 100, n = 1,A = and different r: r = 1.0 (solid line), 
r = 1.498 (dashed line), r = 2.0 (dashed-dotted line). The system was initially in the 
state l^i). 

occurs has a simple interpretation in terms of the covariant matrix elements. Consider 
the threshold in the long time limit in which we may consider the evolution under the 
Markov approximation, but retaining the non-RWA terms. Under this approximation, 
we can put j(t) — > 70 which then allows us to obtain a simple analytical solution for 
the threshold condition for entanglement. 

It is easy to show that the threshold for two mode entanglement occurs at 

V u (t)V 44 (t) = 1/4, (50) 

so that the two modes are entangled when Vu{t)V^{t) < 1/4, otherwise are separable. 
Note that the covariance matrix element Vn(t) is associated with the relaxation free 
modes whereas the element Vu(t) is associated with the mode that is coupled to the 
reservoir and thus undergoes the damping process. Under the Markov approximation, 
we find from equations (I20l - (l26l) that in the long time limit of t ^> 7" 1 , the element 
V44(t) reaches the stationary value equal to the level of the thermal fluctuations 

V M (t) -+2N + 1, (51) 

whereas Vii(t) retains its time dependent behavior which depends on the initial values 

V u (t) = V n (0)cos^ Ft + ^(^d) 2 . (52) 
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We point out that the dependence on the initial values of the long time behavior of 
Vu(t) is due to the presence of the constant of motion V^. 

Averaging equation ( 1521) over a long period of oscillations, the threshold 
condition (1501) simplifies to 

2V n {0) (27V+1) = 1. (53) 

We see that the threshold behavior of entanglement depends on the initial value of the 
covariance matrix element Vii(O). In other words, the entanglement behavior can be 
controlled by the suitable choosing of the initial state. For example, with the initial 
state ( |32|) . we find from equations ( 1531) and ( 136]) that continuous entanglement occurs 
for the degree of squeezing 

r= ~ln(2iV + l) . (54) 

With the parameter value fc^T = 10HQ, we find that the threshold value for r equals 
to 1.498 that is the same found numerically in figure [H We should point out here 
that the same condition for the threshold value of r has been found under the RWA 
approximation [17J. Thus, we may conclude that the threshold value for continuous 
entanglement is not sensitive to the RWA approximation. 

We now proceed to discuss the dependence of the long time entanglement on the 
relaxation rate 70. An example of this feature is shown in figure |2j It is interesting to 
note that under the relaxation the entanglement oscillates in time and the amplitude 
of the oscillations increases with increasing 70 leading to a better entanglement when 
the oscillators are strongly damped. It is a surprising result as one could expect 
that entanglement should decrease with increasing 70 . Again, a straightforward 
interpretation of this effect can be gained from a qualitative inspection of the properties 
of the transformed covariance matrix. 

It is easy to see from equations f )24|) and (125]) that in the limit of vanishing damping, 
j(t) — > and A = 0, the matrix A 2 (i) reduces to Ai(0). One could argue that in 
this limit the covariance matrix elements determined by the matrix A 2 (t) coincidence 
with the elements determined by the matrix Ai(0). Of course, their time evolution is 
determined by the same equations, but there is a subtle difference in their initial values. 
For example, the initial values of the elements whose evolution is determined by the 
matrix Ai(0) are 

*S(0) = \ {MVL\^ ± le 2 "), (55) 

whereas that one determined by the matrix A 2 are 

VS(0) = \ (M^ r ± ie" 5 ") . (56) 

The initial elements are significantly different that what appears as a squeezed 
component in V^(0), the counterpart in V^(0) appears as an anti-squeezed component. 
This is a crucial difference that has a significant effect on the evolution of an 
entanglement. These two contributions cancel each other that results in no oscillations 
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Figure 2. Time evolution of the negativity r;_ and combined variance ((AXi) 2 ) + 
((AY3) 2 ) for A = 100, n = 1, A = 0, r = 1.6 and different values of the relaxation rate 
7o : 7o = 0.05 (solid line), 70 = 1.0 (dashed line), 70 = 5.0 (dashed-dotted line). The 
system was initially in the state \ipi). 



in the entanglement evolution when 70 <C 1. On the other hand, for a large 70 the 
covariance matrix elements determined by A2(t) are rapidly damped to their stationary 
values leaving the elements determined by Ai(0) continuously oscillating in time. These 
oscillations lead to the continuous oscillation of the entanglement seen in figure |2j 

One can interpret these results in terms of collective symmetric and antisymmetric 
states of an iV-atom Dicke model [251 I2S1 |27]. The symmetric and antisymmetric 
states correspond to the atomic dipole moments oscillating in-phase and out-of-phase, 
respectively. The most interesting is that in the case of the atoms coupled to a common 
reservoir, the antisymmetric states do not decay, whereas the symmetric states decay 
with an enhanced rate Nj, where 7 is the single atom decay rate. Hence, in the 
absence of the damping, oscillations induced by the symmetric and antisymmetric states 
cancel each other as they occur with opposite phases. When damping is included, the 
oscillations induced by the symmetric states are damped in time whereas the oscillations 
induced by the antisymmetric states remain unaffected. The oscillations induced by the 
symmetric states die out on the time scale of t ~ 1/ (N'j) leaving the oscillations induced 
by the antisymmetric states unaffected. 

Figure [3] shows the evolution of entanglement and squeezing when the oscillators 
are coupled to each other. In this case there is no continuous stationary entanglement. 
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Figure 3. Time evolution of the negativity r;_ and the combined variance ((AXi) 2 ) + 
((Af 3 ) 2 ) for 7o = 0.05, A = 100, n = 1,A = 0.8 and different r: r = 1.0 (solid line), 
r = 1.498 (dashed line), r = 2.0 (dashed-dotted line). The system was initially in the 
state l^i). 

Thus, the interaction between the oscillators has a destructive effect on the stationary 
entanglement. However, for a large squeezing, entanglement re-appears in some discrete 
periods of time, exhibiting periodic sudden death and revival of entanglement. In 
other words, the threshold behavior of entanglement is a periodic function of time. 
As before, this feature has a simple interpretation in terms of the covariance matrix 
elements. According to equation ( 1501) . for a given temperature the threshold value for 
entanglement depends on the covariance matrix element Vn(t) which, on the other hand, 
depends on A through the frequency parameter Qp- We see from equation (jHJ) that Qf 
decreases with increasing A. Thus, according to equation ( )52|) for interacting oscillators 
the matrix element Vn(t) oscillates slowly in time. The averaging over the oscillations is 
not justified and thus the threshold condition for entanglement is the oscillating function 
of time even in a long time regime. 

Finally, in figure H] we illustrate the evolution of entanglement for two different cases 
of the initial asymmetric state j^)- As we have shown in section [3j more constants of 
motion are then involved than in the symmetric case which, on the other hand, may 
lead to a better stationary entanglement. In the first case, we plot the negativity r]2 
which describes entanglement between the mode 2 and the pair 1 -h- 3. We see from 
figure H(a) that the stationary entanglement appears only when r < r s . Otherwise, the 
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Figure 4. Time evolution of the negativity (a) 772 and (b) 773 for the initial asymmetric 
state ^2) with 70 = 0.05, A = 100, A = 0,r s = 1.489 and different r : r = 1.0 (solid 
line), r = 1.489 (dashed line), r = 2.0 (dashed-dotted line). 

initial entanglement rapidly decays to zero and disappears after a finite time. Again, this 
feature can be easily explained in terms of the transformed oscillators. When r < r s , 
the pair of modes 1 and 2 that is decoupled from the environment is more strongly 
correlated than the pairs 1 -H- 3 and 2 <H- 3, which involve the mode coupled to the 
environment. This preserves the entanglement in the system. In the opposite case of 
r > r s , a large entanglement is initially encoded into the pairs that are damped due 
to the coupling to the environment. This results in the loss of the correlations and 
entanglement. Quite different properties exhibits entanglement between the mode 3, 
which is coupled to the environment, and the remaining pair 1 -h- 2. In this case, 
illustrated in figure H(b) there is no stationary entanglement. This can be interpreted 
as the result of the coupling of the mode 3 to the reservoir that leads to the continuous 
dissipation of the initial correlations r$. 

6. Conclusion 

We have analyzed dynamics of a set of iV harmonic oscillators coupled to a non- 
Markovian reservoir in terms of the covariance matrix. By performing a suitable 
transformation of the position and momentum operators of the system oscillators, we 
have shown that the set of coupled differential equations for the covariance matrix 
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elements splits into decoupled subsets of smaller sizes involving only three and four 
equations. In other words, our analysis clearly show that the dynamics of iV oscillators 
can be completely determined by properties of 4 x 4 and 3x3 matrices. The approach 
proposed here could be particularly useful in applications to macroscopic systems 
composed of a large number of oscillators for which numerical analysis are technically 
too complicated or impossible to perform. 

The approach has been applied to the case of three coupled harmonic oscillators 
interacting with a non-Markovian reservoir. A general feature of the entanglement 
evolution is that it exhibits two characteristic time scales, a shot time regime where an 
initial entanglement is rapidly damped and a long time regime where the entanglement 
undergoes continuous undamped oscillations. Depending on the initial amount of 
entanglement encoded into the system, it can be preserved for all times or may 
periodically disappear and reappear that the entanglement may undergo the sudden 
death and revival phenomena. We have also found that in contrast to what one could 
expect, a stronger damping of the oscillators leads to a better stationary entanglement 
than in the case of a weak damping. Finally, we point out that the three-mode 
entanglement can be observed experimentally, simply by detecting quadrature squeezing 
of the field modes. 
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Appendix A. Initial values of the covariance matrix 



In this appendix, we list the non-zero elements of the initial covariance matrix for the 
case of the asymmetric initial state (1341) . The diagonal elements are of the form 

e -2r s 

V n (0) = -— [9 + e 3rs (3 cosh f - q sinh f )] , 

Q 2r s 

V 22 (0) = — [9 + e~ 3rs (3 cosh f + q sinh f )] , 

e -2r a 



V 33 {0) = —— [l + e 3ra (3 cosh f-q sinh f )] 
8 

e 2r s 

y 44 (0) = — - [l + e" 3rs (3 cosh f + q sinh f )] 



e s 

^55(0) = — (3 coshf + gsinhf) , 
6 

e~ rs 

^66(0) = — — (3 cosh r — gsinhr) , (A.l) 
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whereas the off-diagonal terms are 



_ e -2r s 



Vi3(0) = -j=- [3 — e 3rs (3 cosh r — gsinhr)] , 

8v 3 

2r 

V 24 (0) = [3 - e" 3rs (3 cosh f + q sinh f )] , 

8v 3 

y 35 (o) = Vsmo) = 

v a ( ) = ^(0) = e ' r ' (r °- r - )Bi " hf , 

yor 



where r = a/ 8r% + r 2 s and q = (8r + r s )/r. 



